function  polynom_coeff  = compute_value_function(Y,X)

%Rescaling inputs and output for a more accurate matrix inversion. 
%The estimated polynomial coefficients are then rescaled back;
            
State_var = [ones(size(X(:,1))), log(X(:,1)), log(X(:,2)), log(X(:,1)).*log(X(:,2))];
temp_max = max(State_var,[],1);
X_var_max = repmat(temp_max,[size(State_var,1),1]);
X_var_rescaled = State_var./X_var_max;
Y_max = max(Y);            
Y_rescaled = Y./Y_max;
temp_coeff = (X_var_rescaled'*X_var_rescaled)^(-1)*(X_var_rescaled'*Y_rescaled);
coeff = temp_coeff.*( Y_max./temp_max )';

polynom_coeff = coeff;